#pragma once
#include <Common/headers.h>
#include "pt2curve_projection.h"

namespace shermit {
namespace fitting {
namespace curve {

	/* Cuda Solver of PDM */
	template<typename Float_t>
	class CudaPDMSolverLU {
	public:
		using mat_t = thrust::device_vector<Float_t>;

		CudaPDMSolverLU(
			const std::vector<nv::vec2<Float_t>>& data_pts,
			const std::vector<nv::vec2<Float_t>>& control_pts,
			int max_iter_solver = 200,
			int max_iter_projection = 20, int spr = 256,
			Float_t pt_half_length = 0.05
		)  :max_iter_solver(max_iter_solver),
			max_iter_projection(max_iter_projection),
			spr(spr), projector(data_pts, pt_half_length)
		{
			curve.setControlPts(control_pts.data(), control_pts.size(), DeviceType::CPU);
			nd = data_pts.size();  // number of data points
			nc = curve.nControlPts();   // number of control points
		}

		~CudaPDMSolverLU() {
			curve.deallocate();
		}

		virtual void solve(Float_t error_target) = 0;

		const CubicBSplineCurve<Float_t, DeviceType::GPU> getDeviceCurve() const {
			return curve;
		}

		CubicBSplineCurve<Float_t, DeviceType::CPU> getHostCurve() const {
			auto hcurve = CubicBSplineCurve<Float_t, DeviceType::CPU>();
			curve.copyToHost(hcurve);
			return hcurve;
		}
	protected:
		thrust::device_vector<int> ipiv;
		mat_t pt_soa;
		mat_t workspace;
		mat_t A;
		mat_t tau;    // used for QR factorization
		int nd, nc;

		Pt2CurveProjection<Float_t, DeviceType::GPU> projector;
		CubicBSplineCurve<Float_t, DeviceType::GPU> curve;
		int max_iter_solver;
		int max_iter_projection;
		int spr;
	};
	/*
	class CudaPDMSolverLU32 :public CudaPDMSolverLU<float> {
	public:
		//using CudaPDMSolverLU<float>::CudaPDMSolverLU<float>;
		CudaPDMSolverLU32(
			const std::vector<nv::vec2<float>>& data_pts,
			const std::vector<nv::vec2<float>>& control_pts,
			int max_iter_solver = 200,
			int max_iter_projection = 20, int spr = 256,
			float pt_half_length = 0.05
		):CudaPDMSolverLU<float>(data_pts,control_pts,max_iter_solver,max_iter_projection,spr,pt_half_length)
		{}
		void solve(float error_target = 0.1);
	};
	*/
	class CudaPDMSolverLU64 :public CudaPDMSolverLU<double> {
	public:
		//using CudaPDMSolverLU<double>::CudaPDMSolverLU<double>;
		CudaPDMSolverLU64(
			const std::vector<nv::vec2<double>>& data_pts,
			const std::vector<nv::vec2<double>>& control_pts,
			int max_iter_solver = 200,
			int max_iter_projection = 20, int spr = 256,
			float pt_half_length = 0.05
		) :CudaPDMSolverLU<double>(data_pts, control_pts, max_iter_solver, max_iter_projection, spr, pt_half_length)
		{}
		void solve(double error_target = 0.1);
	};
}
}
}
